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We explore the prospects for constraining cosmology using gravitational-wave (GW) observations 
of neutron star binaries by the proposed Einstein Telescope (ET), exploiting the narrowness of the 
neutron star mass function. This builds on our previous work in the context of advanced-era GW 
detectors. Double neutron-star (DNS) binaries are expected to be one of the first sources detected 
after "first-light" of Advanced LIGO. DNS systems are expected to be detected at a rate of a few 
tens per year in the advanced era but the proposed Einstein Telescope (ET) could catalog tens, if 
not hundreds, of thousands per year. Combining the measured source redshift distributions with 
GW-network distance determinations will permit not only the precision measurement of background 
cosmological parameters, but will provide an insight into the astrophysical properties of these DNS 
systems. Of particular interest will be to probe the distribution of delay times between DNS- 
binary creation and subsequent merger, as well as the evolution of the star-formation rate density 
within ET's detection horizon. Keeping Ho, Slm.o and JIa.o fixed and investigating the precision 
with which the dark-energy equation-of-state parameters could be recovered, we found that with 
10^ detected DNS binaries we could constrain these parameters to an accuracy similar to forecasted 
constraints from future CMB-|-BAO-|-SNIa measurements. Furthermore, modeling the merger delay- 
time distribution as a power-law (cx t") and the star- formation rate (SFR) density as a parametrized 
version of the Porciani and Madau SF2 model, we find that the associated astrophysical parameters 
are constrained to within ~ 10%. All parameter precisions scaled as l/\/iV, where TV is the number 
of cataloged detections. We also investigated how parameter precisions varied with the intrinsic 
underlying properties of the Universe and with the distance reach of the network (which is affected, 
for instance, by the low-frequency cutoff of the detector). We also consider various sources of distance 
measurement errors in the third-generation era, and how these can be folded into the analysis. 

PACS numbers: 98.80.Es, 04.30. Tv, 04.80.Nn, 95.85.Sz 



I. INTRODUCTION 

The era of advanced gravitational-wave (GW) detec- 
tors is approaching quickly. The previous decade has 
seen significant improvements in the sensitivity of GW 
interferometers, leading to the construction and opera- 
tion of two Laser Interferometer Gravitational- wave Ob- 
servatory Q^IGO) [l| detectors in the USA, GEO-600 in 
Germany Virgo in Italy Q and TAMA-300 in Japan 
The latter detector was designed as a testbed to 
develop new technologies for the proposed underground, 
cryogenically cooled KAGRA (formerly LCGT [5]) detec- 
tor ^. The LIGO, Virgo and GEO-600 detectors have 
conducted joint searches since 2007. 

The most promising source for the first detection 
of gravitational waves is the inspiral and merger of a 
compact-object binary consisting of neutron stars (NSs) 
and/or black holes J]- The first joint search for compact 
binary coalescence signals during the LIGO S5 science 
run and the Virgo VSRl data did not result in direct de- 
tections [§[ , nor did the "enhanced" detector search dur- 
ing the LIGO S6 science run and the Virgo VSR2-I-3 data 
@. Furthermore, the upper limits placed on compact- 
binary coalescence rates from the latter search remain 
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two to 3 orders of magnitude above existing astrophysi- 
cally predicted rates. However, the LIGO detectors are 
currently being upgraded to their "advanced" configu- 
ration [10[, due for completion in ~ 2015, for which 
the horizon distance for NS-NS inspiral detection will 
be boosted to ~ 450 Mpc, giving an almost thousandfold 
gain in volume sensitivity of the detectors. The advanced 
detectors are expected to detect double NS inspirals at 
a rate of ^ 40 yr^^, although this may vary by approxi- 
mately 2 orders of magnitude in either direction 

Complementing AdLIGO will be a global network of 
advanced detectors, including AdVirgo [l^, KAGRA 
and possibly a third LIGO detector in India, LIGO- 
India ^i^. There are currently no prospects for a South- 
ern Hemisphere GW interferometer operating in the ad- 
vanced era. A global network comprising these detectors 
will help turn the field from the search for the first de- 
tection, into a precise astronomical tool. 

The GWs emitted by a compact binary system directly 
encode the redshifted masses and luminosity distance of 
the system. Double NS (DNS) binary systems are com- 
monly referred to as self- calibrating standard sirens be- 
cause their distance is directly encoded in the waveform, 
and a means of determining their redshift would mean 
we could probe the cosmic distance ladder and extract 
cosmological parameters (l3 - [T7| . While the phase evo- 
lution of the strain signal in a single interferometer can 
give precise constraints on the redshifted mass of the sys- 
tem, we require a global network of detectors to constrain 
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the sky location, orbital inclination and GW polarization 
so that we can break the degeneracies of these angular 
factors with the luminosity distance. 

Unfortunately the redshift and intrinsic mass of the 
systems enter the waveform only in a combination as the 
redshifted mass parameter; hence previous techniques for 
performing GW cosmology using these standard sirens 
has relied on the association of the GW source with an 
electromagnetic counterpart, which can provide an inde- 
pendent redshift measurement [IStSII]. In our previous 
paper, Taylor et al. [l^], we studied a technique for prob- 
ing the Hubble constant and NS mass-distribution pa- 
rameters using only the GWs detected by an advanced- 
era network. This work relied on the assumption, sup- 
ported by observations, that the NS mass-distribution is 
sufficiently narrow, which means that we already have 
a good idea of the intrinsic masses in the system and 
the measured redshifted mass parameter then provides a 
narrow distribution of possible redshifts for each observed 
source. Combining these redshift distributions with the 
network-measured luminosity distance for a catalog of 
observed DNS systems provides constraints on the un- 
derlying cosmological parameters, as well as the astro- 
physical distribution of these systems. This technique 
(first considered by Markovic in Ref. [2^ and extended 
in Refs. [U,!^^) relies on measurements of the redshifted 
chirp mass (expected to be the best-determined param- 
eter) and luminosity distance for a catalog of detected 
systems. 

In our previous analysis, the cosmological parame- 
ters that we could constrain were restricted by the sub- 
Gpc reach of an advanced-era network. We now ex- 
tend this technique to a third-generation network, which 
could have a reach out to tens of Gpcs. Proposed third- 
generation detectors aim for a broadband factor of 10 sen- 
sitivity improvement with respect to advanced detectors, 
and to increase sensitivity in the range '-^ 1 — 10 Hz, com- 
pared to the 10 — 20 Hz lower frequency cutoff of ad- 
vanced detectors. As a prototypical third-generation de- 
tector we use the Einstein Telescope, consisting of three 
overlapping 10 km arm-length interferometers in a tri- 
angular configuration [26l428l |. Each interferometer may 
actually be two detectors: a cryogenically cooled, under- 
ground detector with good low-frequency sensitivity, and 
a high laser power detector with good high-frequency sen- 
sitivity [2^. Keeping Ho, ilm.o and JIa^q fixed, we find 
that the sensitivity provided by such a network will be 
large enough to constrain the dark-energy equation-of- 
state parameters and NS mass-distribution parameters, 
as well as the astrophysical distribution of the systems. 
The latter will inform us about the average time delay 
between the formation of these compact-binary systems 
and their merger, as well as the shape of the underlying 
star-formation-rate density. 

Third-generation detectors are unlikely to be online 
before the mid-2020s, but, if realized, the ambitious and 
novel design for the Einstein Telescope will have far- 
reaching scientific advantages. Such a detector could de- 



tect as many as hundreds of thousands of DNS inspi- 
rals per year, which, along with the distance reach of 
the detectors, will permit precision GW astronomy. In 
this paper, we will not consider other methods that have 
been proposed for using GW observations as cosmologi- 
cal probes. In particular, we do not consider association 
of GW detections with an electromagnetic (EM) coun- 
terpart, which has been studied in Refs. p^. Il7|. nor do 
we consider tidal-coupling corrections to the phase evo- 
lution of the strain signal [30| . The latter method is also 
a GW-only technique with significant potential, in that 
these phase-evolution corrections break the mass-redshift 
degeneracy at 5PN order and, assuming the NS equation- 
of-state is well known, will permit the distance-redshift 
relation to be probed. It may also be possible to ap- 
ply the method used by Ref. [31| , which was investigated 
in the context of future space-based detectors, to third- 
generation ground-based detectors. Their method relies 
on the measurement of cosmologically-induced shifts in 
the GW-phase at 4PN order. 

This paper is laid out as follows. Section HIl describes 
the characteristics of the Einstein Telescope, as well as 
possible third-generation networks and detection criteria. 
In Sec. HIT] we discuss aspects of DNS systems, including 
the mass distribution of the constituent NSs, and the red- 
shift distribution of DNS mergers. Section IIVI describes 
the effect of the dark-energy equation-of-state parameter 
on cataloged luminosity distances, while Sec.lVldescribes 
how we construct and analyze catalogs of detected DNS- 
system inspirals. Results are shown in Sec. IVIl followed 
by our conclusions in Section [VIII 



II. DETECTOR CHARACTERISTICS AND 
NETWORKS 

A. The Einstein Telescope 

The Einstein Telescope (ET) is a proposed third- 
generation ground-based interferometer. A recent design 
study has been carried out within the European Commis- 
sion's FP7 framework [35| to evaluate the science case 
for such a detector, and to consider the technological 
advances required for the science goals to be achieved. 
Through this three-year design study, some favored de- 
signs and configurations have emerged. 

The aim for third-generation ground-based detector 
designs is to achieve a broadband factor of 10 sensitivity 
improvement with respect to advanced detectors, and to 
push the sensitivity down into the ^ 1 — 10 Hz range. 
Early designs examined the prospects for pushing con- 
ventional techniques used in advanced detectors to their 
limits to construct a third-generation interferometer [SSj . 
This gave the ET-B noise curve in Fig. [T] Beyond the 
extension of the arm-length to 10 km, several techniques 
were proposed to suppress high- and low-frequency noise, 
including siting the detector underground. 

Crucially, the techniques used to suppress high- 
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FIG. 1: Comparison of the different noise curves for AdLIGO 
(high-power-zero-detuning) ^32i], the initial Einstein Telescope 
noise curve based on conventional techniques, ET-B |33| , the 
initial "xylophone" noise curve, ET-C [34l , and the improved, 
more realistic xylophone noise curve, ET-D |29|| . These noise 
curves are for one 10km right-angled interferometer. 



frequency noise are not necessarily compatible with the 
suppression of low-frequency noise. Increasing the laser 
power will reduce the photon shot noise which dominates 
the high-frequency range, but this worsens the thermal 
noise which dominates at low frequencies. The "xylo- 
phone" design was proposed to address this issue. In- 
stead of having a single broadband instrument, the xy- 
lophone design comprises a high-power, high-frequency 
interferometer (ET-HF) and a cryogenic low-power, low- 
frequency interferometer (ET-LF) \M\. ET-LF would be 
an underground instrument, and limited at low frequen- 
cies by gravity-gradient noise, while ET-HF would be 
colocated and co-oriented with ET-LF but surface-sited. 
ET-HF would employ high-power lasers to suppress high- 
frequency photon shot noise. 

The initial xylophone design gave the ET-C curve [s^l 
in Fig. [T] which was refined to give the ET-D xylophone 
design [29|]. We will use the ET-D noise curve in the 
ensuing analysis. 

Current and advanced era ground-based interferom- 
eters are right-angled interferometers, since an arm 
opening-angle of 90° maximizes their sensitivity. How- 
ever, if both GW polarization states are to be measured 
at a single site, then two or more colocated nonaligned 
interferometers are required. Furthermore, at least three 
colocated interferometers are required to construct a null- 
stream, i.e., a sum of individual interferometer responses 
that is insensitive to GWs and can be used to identify 
noise transients in the data stream. 

Taking these goals into account, the design requiring 
the shortest total length of tunnels is a triangular con- 
figuration with three identical interferometers positioned 
at each vertex of the triangle, an arm-openin g angl e of 
60° and rotated relative to each other by 120° [26l - l28j . A 



triangular configuration also provides redundancy, since 
polarization constraints are still possible with only two 
vertices operational. 

In the following we consider three ET-D interferome- 
ters in the triangular configuration, which we denote as a 
"single FT" . More than one ET would be very optimistic, 
so we complement our single ET with a network of indi- 
vidual third-generation right-angled interferometers (also 
with ET-D sensitivity) to permit source distance deter- 
mination. While different locations are being mooted, we 
choose the Virgo location as the reference ET site |3a |. 



B. Signal-to-noise ratio 

The optimal matched filtering signal-to-noise ratio of 
a GW detection is given by 



Popt 



^ S,{f) 
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(1) 



where Sh{f) denotes the detector's noise power spectral 
density. In the quadrupolar approximation, the Fourier 
transform of the signal amplitude of GWs from an inspi- 
raling binary system takes the following form [s^, HI] : 
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The function Q is defined by 

e = 2[Fl{l + COS^ t)2 + 4i^2 ^Qg2 ^]l/2^ 

where < < 4, and 

= -(1 + cos^ 6) cos 20 cos 2'(/; — cos6'sin2(/)sin2'(/;. 

Fx = -(1 + cos^ 6') cos20sin2V' + cos6'sin20cos2'(/; 

^ (4) 

are the interferometer's strain responses to the different 
GW polarizations. 

Following Ref. [1^ we can write the matched filtering 
signal-to-noise ratio in a single detector as 
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where Aiz = (1 + is the redshifted chirp mass. 
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TABLE I: The characteristic distance sensitivities [as given by 
evaluating Eq. ^] of some advanced-detector configurations and 
various design studies for the third-generation Einstein Telescope. 



Detector 



ro / Mpc 



AdLIGO" 
AdLIGO' 
AdLIGO" 
AdVirgo'' 
ET" 
ET-B^ 
ET-C» 
ET-D'' 



80 
110 
119 

85 

1527 
1587 
1918 
1591 



"Low-power zero detuning 3^ 
'High-power zero detuning | 3^ 
'^NS-NS optimised ^ 
■^Ref. [H 

'^Polynomial noise-curve approximation [4^ ] 
■'^Conventional technology 35j 

53'''^-generation technology, xylophone configuration '3E?| 
''3'^'^-generation technology, xylophone configuration (updated 
and more realistic) [35l | 



and 2/i„ax is the wave frequency at which the inspiral 
detection template ends [39|. The intrinsic chirp mass, 
M , is given in terms of the component masses by, 
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M = 



77117712 



(mi + 7)72) 



{nil + rri2)- 
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C. Network antenna patterns 

The angular dependence of the SNR is encapsulated by 
the variable 0. The sky location and binary orientation 
can be deduced from the network analysis, however in 
our analysis we will use only Dl and A4z measurements. 
We calculate the probability density function for Q |25j | 
numerically using Eq. ^ by choosing cos 9, (fi/n, cost 
and ip/iT to be uncorrelated and distributed uniformly 
over the range [—1,1]. 

It is unlikely that more than one ET will be con- 
structed. A more likely network configuration will involve 
a single ET with single third-generation right-angled de- 
tectors at other sites. In the interest of verisimilitude 
we take into account possible detector locations for such 
a third-generation network. Table [ll] contains the loca- 
tions and arm-bisector orientations of various current and 
planned detectors. 

To write down the antenna pattern function as a func- 
tion of the detect or p osition,^ we use the notation and 
formalism of Ref . [43| . 

For a GW source at coordinates {6, (f)) on the sky. with 
polarization angle -0 and a detector with opening angle t] 
at latitude /? and longitude A and such that the bisector 
of its arms points at an angle x^ counterclockwise from 
East, the antenna pattern functions are 




= svarj 



cos {^Tp) sin (2?/;) 
— sin {2i}}) cos (2-0 




(9) 



where, 



The phase evolution of the strain signal in a single in- 
terferometer can constrain A^^ to subpercent precision 
(ioj . In order to measure the luminosity distance, Dl, 
we require a network of separated detectors to break the 
waveform degeneracy between Dl and Q [see Eq. ([2])]. 
Distance measurement errors in a third-generation net- 
work will be discussed in Sec. IVIDI 

The signal-to-noise (SNR) of a detected system will 
vary between the individual network sites, as a result of 
the different ShifYs and angular dependencies. However, 
following Refs. [23, [1^, we assume the network SNR of 
a detected system is given by the quadrature summation 
of the individual interferometer SNRs, 



Pnct 




(ro,fcefcCfe(/max))'. (8) 



where ro.fe, Cfc(/max) and Qk encode the distance, fre- 
quency and angular sensitivities of the fc*^ detector. A 
comparison of the characteristic distance sensitivities of 
some second- and third-generation detectors is shown in 
Table H 



a=— sin (2x) [3 - cos (2/3)] [3 - cos (261)] cos [2(0 + A)] 
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cos (2x) sin/3[3 - cos (261)] sin [2(0 + A)] 



-I- - sin (2x) sin (2/3) sin (29) cos (0 + A) 



cos (2x) cos /3 sin (29) sin (0 + A) 
sin (2x) cos^ /3 sin^ 9, 



b — cos (2x) sin /3 cos 6 cos [2(0 + A)] 
1 



■ sin (2x) [3 - cos (2;3)] cos 9 sin [2(0 + A)] 
COS (2%) COS /3 sin 9 cos (0 + A) 
i sin (2x) sin (2/3) sin 9 sin (0 + A) . 



(10) 



As a reference, we use a network comprising three 60° 
ET-D sensitivity interferometers at the Virgo location 
(a single ET), plus right-angled interferometers at the 



^ We do not consider modulation of the antenna patterns due to 
the rotation of the Earth. We justify this in Sec. IVI Fl 
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TABLE II: A reproduction of the GW-interferometer geographical locations, and arm-bisector orientations from Schutz [43| . We 
include updated IndlGO information 44]. 



Detector 



Label 



Longitude 



Latitude 



Orientation 



LIGO Livingston, LA, USA 
LIGO Hanford, WA, USA 
Virgo, Italy 

KAGRA (formerly LCGT), Japan 
LIGO-India, India 



L 
H 
V 
J 
I 



90°46'27.3" W 
119°24'27.6" W 
10° 30' 16" E 
137°10'48" E 
76° 26' E 



30°33'46.4" N 
46°27'18.5" N 
43°37'53" N 
36°15'00" N 
14° 14' N 



208.0° (WSW) 
279.0° (NW) 
333.5° (NNW) 
20.0° (WNW) 
45.0° (NE) 



LIGO-Livingston and LIGO-India locations. The char- 
acteristic distance reach of all of the interferometers in 
the network is taken as 1591 Mpc, corresponding to ET-D 
sensitivity [1^. This is the sensitivity of a 10 km right- 
angle interferometer. We account for the different detec- 
tor arm-opening angles in the antenna pattern functions. 

The network SNR given by Eq. ^ also depends on 
C(/max), which describes the overlap of the signal power 
with the detector bandwidth [2^. The frequency at the 
end of the inspiral (taken to correspond to the innermost 
stable circular orbit) is at 



785 Hz /2.8M, 



1 + z 



M 



(11) 



where AI is the total mass of the binary system |37| . 
The maximum binary-system mass could conceivably be 
~ 4.2Mq.^ The ET horizon distance for a system with 
a total mass of ~ 4Mq is ~ 25 Gpc [la]- In the 
ACDM cosmology this corresponds to a redshift of ~ 2.9, 
and from Eq. (in]) this gives /max ^ 134 Hz. Given 
the ET-D noise curve y^C(/max = 134Hz) > 0.98. 
E xtending the red shift reach out to z ^ 5 still gives 
VC(/max = 87Hz) > 0.96. Thus, we feel justified in 
adopting C(/inax) — 1 for all interferometers in the en- 
suing analysis. 

Using these expressions we were able to numerically 
estimate the probability distribution for the effective 0, 



e 



off 



(12) 



where the sum is over all detectors in the network. We 
use this 8off distribution to choose SNRs for each source 
in the catalog via Eq. ([SJ and then impose a detection 
criterion. As a reference, we adopt the detection criterion 
that the network SNR must be greater than 8. 



III. DNS SYSTEMS 
A. Neutron-star mass distribution 

For a full discussion of our assumptions and modeling 
details of the NS mass distribution in DNS systems, see 
our previous work (Ref. j2^ and references therein) . We 
provide here a brief summary of the main assumptions 
pertinent to the present study. 

To lowest order, the GW signal depends on the two 
neutron-star masses through the chirp mass, Ai. We 
assume that the distribution of individual neutron-star 
masses is normal, as suggested by analysis of Galactic 
DNS systems '46^, and population synthesis studies 
(see, e.g., Refs. (3^41, 48]). For ctns ^ Mns, this should 
also lead to an approximately normal distribution for the 
chirp mass. 

We use a simple ansatz for the relationship between 
the chirp-mass distribution parameters and the underly- 
ing neutron-star mass distribution. The chirp mass dis- 
tribution is modeled as normal, 

M^N(pl,,cjI), 

with mean and standard deviation 

^Ji, « 2(0.25)3/5^NS, (Tc ~ y2(0.25)3/5aNs, (13) 

where /iNS and ctns are the mean and standard deviation 
of the underlying neutron-star mass distribution. 

A recent study by Ozel et al. [i^ has found that 
DNS data are consistent with both pulsar and compan- 
ion having been drawn from the same underlying distri- 
bution of masses. The literature indicates an underly- 
ing neutron-star mass distribution in DNS systems with 
CTNS < 0.15Mq [Ulii, liil.^ Hence, we anticipate that 
Eq. will be appropriate for generating data sets and 
we use this in the ensuing analysis. The assumption 
throughout is that for the volume of the Universe probed 
by our global network, the neutron-star mass distribution 
does not change. 



^ Both neutron stars in the binary system would need to have 
masses 2cr above the distribution mean at the maximum consid- 
ered fi and cr, where /^ns S [1.0, I.SJMq, aj^g S [0, O.3]M0. 



^ Indeed, Ozel et al. [4d| indicates that the DNS mass distribution 
is peculiar, since it cannot be explained via electron-capture or 
core-collapse supernovae mechanisms; rather, its narrow disper- 
sion may be a result of the evolutionary path of these systems. 



6 



Further population synthesis and observational studies 
in the following decade will help to shed further light on 
the nature of the NS mass distribution. The assumption 
of a unimodal (for DNS systems) Gaussian distribution 
is an approximation, and if future studies show this to 
be inappropriate, then a more suitable ansatz could be 
readily incorporated within the framework described in 
this paper. 



B. Merger-rate density 

In this analysis, we aim to probe not only the back- 
ground cosmology and NS mass-distribution parameters, 
but the astrophysical properties of the binary NS pop- 
ulation. To this end, we now consider the factors con- 
tributing to the coalescence of a binary NS system. 

Following several population synthesis studies (e.g.. 



Refs. ,38, 50]) and Ref. 
per comoving volume as 



51| . we define the merger rate 



n{t) 



dt 



tb)^^itb)dtb, 



(14) 



where A is a mass efficiency, defined as the number of 
coalescing DNS binaries per unit star- forming mass [SS*] . 
dPjn/dt is the probability density distribution of a DNS 
binary merging at a time (t — tb) after formation, and 
dpf,/dt is the star- formation rate (SFR) density at cos- 
mological time tb- 

Star formation, and the efficiency of double compact- 
object formation from the progenitor system, may be 
sensitive to the host-galaxy morphology and environ- 
ment (e.g. metallicity) . Furthermore, the distribution 
of delay times between star formation and the corre- 
sponding DNS-system coalescence may have contribu- 
tions from several different evolutionary paths. However, 
we are interested here only in a third-generation GW- 
interferometer network's ability to constrain various as- 
trophysical and cosmological parameters [l^j- As such 
we consider a single component star-formation distribu- 
tion, delay-time distribution and mass efficiency, defer- 
ring considerations of the other dependencies to a future 
study. We now discuss the various terms in Eq. ([T4|) in 
more detail. 



1. Mass efficiency, A 

We use values for A obtained from the population syn- 
thesis calculations of Ref. [50] . Smoothed histograms of 
the mass efficiency are shown in Fig. 4 of that paper, 
with modes at ^ 10~^Mq^ for DNS systems formed in 
both elliptical and spiral conditions. However the distri- 
bution of A ranges over several orders of magnitude, with 
IO-^A/qI < a < 10-3Mq\ We adopt A = IQ-^M^^ as 
the reference value for our analysis. 



2. Merger-delay distribution, dPm/dt 

Massive stars in high-mass binary systems evolve into 
DNS systems on much shorter timescales than typical 
galaxy evolution or Hubble timescales, such that dPm/dt 
is essentially determined by the initial orbital separation, 
a, of the DNS system [s^l- The evolutionary time delay 
between the progenitor formation and the formation of 
the corresponding DNS system is typically < 50 Myr [5^ , 
and is therefore negligible compare to the gravitational- 
wave inspiral timescale, which scales as Tgr oc a'*. As- 
suming the number of binaries, A'^, born with separation 
a scales as dN/da oc a'^ [s^ , we obtain 



dPm 

dt 



dN dN da 



dTe-r da dT„ 



(15) 



If DNS systems have the same orbital separation distribu- 
tion as normal-abundance main-sequence stars [s^ . [55| , 
then 7 = a = — 1. However, this scaling is not well 
constrained and this is discussed in more detail in Ap- 
pendix lA 11 Instead, we adopt the approach of allowing 
a to be a free parameter that we attempt to fit from 
our observations and ask with what precision this can 
be determined. We use the value a — —1 for our refer- 
ence model, which is justified by current (albeit sparse) 
analysis of Galactic DNS systems ISa-ISa , and various 
population synthesis calculations [SCt l53l [59l - [6^ . For 
normalisation purposes, we assume a minimum delay- 
time of 50 Myr and a maximum delay-time equal to the 
cosmology-dependent age of the Universe; these choices 
are discussed in Appendix lA II 



3. Star-formation rate density, dpt/dt 

The star-formation rate density is also rather uncer- 
tain. The SF2 model of Porciani and Madau [1^ at- 
tempts to factor in the uncertainties in the incomplete- 
ness of data sets and the amount of dust extinction at 
early epochs. The SF2 model has the form 



~dt 



(^) 



0.16 X 



E{z) 



(l-^z)3/2 



exp {(iiz) 
exp (/32z) + 22 



MqMpc-V"\ 



(16) 



where 



E{z) - ^amdl + zf + fifc,o(l + zY + r!A(z) (17) 

and fi\ — fi2 — 3.4. In this model, the SFR density 
remains roughly constant at z > 2, which may be in- 
compatible with recent observations [13, [6^1 ■ This is dis- 
cussed in more detail in Appendix I A 21 To allow for some 
uncertainty, we treat fi\ and fii as free parameters and 
explore how precisely we can measure them. While this 
simple ansatz does not cover all possible forms for the 
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SFR density, using it will provide an indication of what 
GW observations could tell us. The framework is easily 
adaptable to more complex SFR models. 

We must also specify , the lower integration bound in 
Eq. p4)) . which represents the time of the earliest period 
of star formation. The highest redshift objects observed 
are a long gamma-ray burst (GRB) with a photometric 
redshift of ~ 9.4 fOBs] and a candidate z ~ 10 galaxy 
67| . We therefore use z = 10 as the earliest time of star 
formation. Future observations, for instance with the 
James Webb Space Telescope, may be able to probe back 
to the first phases of galaxy formation at z 15 and if 
objects are found at that epoch, this assumption should 
be revised. However, our results are fairly insensitive to 
this choice. 



4- Calculating h(z) 
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FIG. 2: Merger-rate densities computed for the reference 
cosmology (solid line) and for cosmological parameters chosen 
randomly from within the prior range (red crosses). 



Eq. (jl4p can be rewritten as a distribution in redshift 
using dt/dz = -1/((1 + z)E{z)Hq) 



h{z) 



A^— - h) — [tb)—dzb 



'10 

0.16A 



dt 



10 



dP„ 



Hq J ^ dt 
( exp(/3iZb) 



dt dzb 

{t{z) - tb{zb)) 

dzb 



Vexp(/32Zb) + 22; (l + Zb)5/2 



(18) 



Evaluating this expression requires an expensive double 
integral which created a bottleneck in our analysis. How- 
ever, because the priors on the cosmological parameters 
are narrow (see Sec. IV Pp . there is little variation in the 
merger-rate density across this range, as shown in Fig. 
[2j We therefore fixed the cosmological parameters at 
their reference values for the cosmological time calcu- 
lation, which made the merger-rate density calculation 
considerably faster. Although this throws away some of 
the cosmological information, it did not significantly af- 
fect the results and made the analysis more tractable. 



IV. COSMOLOGICAL MODELS 

In our previous analysis [l^l we considered only a flat 
cosmology, but here we allow for curvature and an evolv- 
ing equation of state (EOS) for the dark energy. From 
the cosmological field equations we have 



p + 3 



(19) 



where p and p are the density and pressure of a cosmo- 
logical fluid respectively, while a is the scale factor of the 
universe. Derivatives are with respect to physical time. 



For a perfect fluid {p = wpc^ 
rameter), this reduces to 



where w is the EOS pa- 



^j=-3(l + .;)(^^). (20) 



Hence, 

p(a) = p(atoday - 1) X e-3/i"(i+"')('^-7a'). (21) 

The last decade has seen many proposals for w, with 
different physical motivation. One approach attempts to 
explain dark energy as a minimally coupled scalar field 
("quintessence") slowly rolling down its potential such 
that it can exert negative pressure. While it is possible to 
have a constant EOS in this formalism, this requirement 
places severe constraints on the potential and so it is 
natural to expect a time- varying EOS |68| . 

The simplest approximation is to assume a linear 
model (w(z) = wq + wiz), but this is only appropriate 
for local studies due to the divergence at high redshift. 
The Shafieloo-Sahni-Starobinsky ansatz [6^ models the 
EOS evolution as a "tanh" form that ensures w — — \ at 
early times and w — ^ at low z. This ansatz prevents the 
crossing of the "phantom divide" at w = —1, desirable 
since phantom fluids can not be explained by a mini- 
mally coupled scalar field The ansatz we adopt in 
this work is the Chevallier-Polarski-Linder ansatz [H, [t^ 



w(a) = ii;o-|-Wa(l — a), w{z) = wo+Wa 



l + z 



(22) 



This ansatz was adopted by the Dark Energy Task Force 
[7l| . and has several desirable features. It depends on 
only two free parameters, it reduces to the linear model 
at low z, and it is well-behaved at high redshift, tending 
to wo + Wa- Using this EOS, 



17a(z) = f^A^o X (1 + z) 



(23) 
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For different global geometries of the Universe the lu- 
minosity distance, D^, is given by 



where 



T{z\C) 




^k,Q > 0, 
^k,Q = 0, 
^k,0 < 0, 



(24) 

in which Dh is the Hubble scale (c/Hq) and 
C={Ho, fim,o, ^A,o, ^k,o, WQ,Wa} is the set of cosmologi- 
cal parameters describing the large-scale characteristics 
of the universe. 

The comoving radial distance, Dc{z), is given by 



Dc{z) = Dh 



dz' 

wy 



(25) 



where E{z) is given by Eq. (IT71) . The redshift derivative 
of the comoving volume is given generally by 



dz "{i + zYE{zy 



(26) 



V. MAKING & ANALYZING DNS CATALOGS 

We refer the reader to our previous study [22!| for full 
details of our calculation, but we summarise the main 
details here. 

A. Distribution of detectable DNS systems 

The two system properties we will use in our analysis 
are the redshifted chirp mass, Aizi and the luminosity 
distance, D^. We assume that only systems with an SNR 
greater than a given threshold will be detected. We can 
write down the distribution of the number of events per 
unit time in the observer's frame with Ai, z and effective 

e [Ulli, 



d'^N 



dVc h{z) 
dtdQdzdM dz (1 + z) 



V{M)Ve{Q). (27) 



The 1/(1 -|- z) factor accounts for the redshifting of the 
merger rate [SSj . 

Converting this to a distribution in Mz, Dl and p, and 
integrating over p to find the distribution of detectable 
systems (i.e., systems above SNR threshold) gives 



d^N 



dtdDLdMz 



P>Po 



dz dVc n{z) 
dDL dz (1 -I- zy 



X V 



Mz 



D, 



X Ce 



Po Dl 
8 ro 



1.2AfQ 



5/6' 



(28) 



where the form of (dz/dD^) will depend on the curvature 
of the Universe [see Eq. (|24| ]. 

To calculate the number of detected systems (given 
a set of model parameters, 7^) we integrate over this 
distribution, which is equivalent to integrating the dis- 
tribution over redshift and chirp mass, i.e. N/^^ = T x 

dzdAi, where T is the duration of the 



rOO rOO / 

Jo Jo \'<: 



d-'N 



dtdzdAi 

observation run.* 



B. Creating mock catalogs of DNS binary 
inspiraling systems 



The model parameter space we investigate is the 7D 
space of [wo, Wa, /iNSj ctns, /32]- To generate a cat- 
alog of events, we choose a set of reference parameters, 
motivated by previous analysis in the literature. For our 
reference cosmology, we adopt Hq = 70.4 kms~"'^Mpc~^, 
p^j^o = 0.2726, r^A.o = 0.728, wq = -1.0 and Wa = 0.0 
[72] ■ The reference parameters of the neutron-star mass 
distribution are //ns = 1.35Mq and ctns = O.O6M0 [isl ]. 
We have previously discussed the delay-time distribution 
and SFR density in Sec. IIIIBI We adopt a power-law 
merger-delay distribution with reference power-law index 
a = —1.0, and we take the SFR density to be given by 
the SF2 ansatz [63], with /3i = 3.4 and ^2 = 3.4. 

These reference parameters are used to calculate an 
expected number of events^, and the number of observed 
events. No, is drawn from a Poisson distribution (assum- 
ing each binary system is independent of all others) with 
that mean. Monte-Carlo acceptance/rejection sampling 
is used to draw random redshifts and chirp masses from 
the distribution in Eq. 1^7^ for each event. The Dl and 
Mz are then computed from the sampled M and z. 

For these reference parameters, which give a local 
merger-rate density of ~ 2 x 10~^ Mpc~^yr~^, and 
assuming detected systems must have a network SNR 
greater than 8, we find that the expected number of de- 
tections in 1 yr is ^ lO'^.^ 



We found that, for the purposes of the calculation, assuming the 
NS mass distribution was a (5-function, centred at the mean given 
by the trial parameters, allowed at least a tenfold speedup in the 
calculation. See Appendix [B] for further details. 
The observation time, T, is assumed to be 1 yr, and the mass 
efficiency is assumed to be A = 10~^Mg ^ (as mentioned earlier). 
This reference local merger-rate density is roughly five times 
smaller than the realistic value quoted in Abadie et al. [ill , but 
20 times larger than the pessimistic estimate. Whilst we could 
scale our merger-rate density calculations to match the realis- 
tic value of 10~® Mpc~''yr~^, our modified likelihood statistic 
makes our analysis insensitive to such scalings. A rescaling to the 
realistic local merger-rate density of Ref. [llj would lead to an 
expected detection rate of approximately half a million sources. 
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C. Likelihood statistic 

In the measurement-parameter space of redshifted 
chirp mass and luminosity distance, the measured num- 
ber of detections in a given bin will be a Poisson random 
variate with a model-dependent mean. If we take the 
continuum limit of this, such that bin sizes are infinitesi- 
mally small and contain either or 1 events, then we can 
formulate the likelihood of a catalog of discrete events, 



/:(5|7t,H)- 



No 

^n^(Ati7^), 



(29) 



where A 



{At, 



with Xi = 



IS the vector of measured 
system properties, with = {Mz,DL)i for system i. 
No is the actual number of detected systems, while Nfj_ is 
the number of DNS inspiral detections predicted by the 
model with parameters 7^. Finally, r(Ai |7^) is the rate of 
events with properties A4z and Dl, evaluated for the i^^ 
detection under model parameters 7^, which is given by 
Eq. (pS)) . The trial cosmological parameters are used to 
calculate a model-dependent redshift from the cataloged 
luminosity distance, and, in turn, this redshift is used 
to infer a model-dependent intrinsic chirp mass from the 
cataloged value of Mz- These values oi Mz, Dl, as well 
as the model-dependent values of z and Ai are inserted 
into Eq. (|28|) to compute the likelihood. 

In our previous study, we employed a modified likeli- 
hood statistic which marginalized over the local merger- 
rate density of DNS systems. This approach reflects our 
current lack of knowledge of this quantity, estimates of 
which vary over several orders of magnitude. We adopt 
the same approach in this analysis, to eliminate the de- 
pendence on poorly known scaling factors. This includes 
the mass efficiency parameter. A, which is a quantity de- 
rived from population synthesis studies and can vary over 
several orders of magnitude. 

The modified likelihood statistic is 



C(i\-ji,H)<xN,, 



No 



(30) 



1=1 



We note that we have not included a prior on the scal- 
ing factors in the above, which is equivalent to using an 
improper flat prior over the range [0, oo]. This reflects 
our current lack of knowledge of the mass efficiency. 



D. Calculating the posterior probability 

We use a weakly informative prior on the model pa- 
rameters, so that it does not prejudice our analysis. As a 
prior on /iNS "we take a normal distribution with parame- 
ters /i = I.35M0, a = O.13M0. This is motivated by the 
posterior predictive density estimate for a neutron star 
in a DNS binary system given in Ref. (45| . 



Given that ET will most likely not be operational until 
the mid-2020s, we must consider what constraints con- 
ventional observational cosmology can put on cosmolog- 
ical parameters. In the recent study by Zhao et al. jl7| . 
the authors investigated how the dark-energy EOS could 
be probed by ET observations of DNS systems to provide 
distance measurements, complemented by electromag- 
netic measurements of the associated short gamma-ray 
burst (sGRB) to provide the redshift. They estimated 
that a combination of the Planck cosmic microwave back- 
ground (CMB) prior, JDEM BAO resuhs, and future 
Type la supernova observations could provide cosmolog- 
ical constraints by the ET era of 

An,nfi = 3.46 X 10"^, Arjfe,o = 5.91 x 10^", 
AHo = 0.336, Awo =0.048, Awa = 0.184. (31) 

Hence, we assume that Hq, ilm.o and flk.o are precisely 
known quantities, with values of 70.4 kms~^Mpc~^, 
0.2726 and —0.0006, respectively. As a prior on wq and 
Wa, we adopt the constraint that w{z) < —1/3 at all 
redshifts. Hence we use uniform priors on the EOS pa- 
rameters with Wq < —1/3 and WQ + Wa < ~l/3 and lower 
bounds set low enough so as not to affect the posterior 
probability distribution. We also adopt uniform priors 
for all other parameters under investigation. 



E. Bayesian analysis and an adaptive Markov chain 
Monte Carlo technique 

Bayes' theorem states that the posterior probability 
distribution of the parameters 7^ describing a hypothesis 
model H, and given data D is given by 



^' ' ' ^ p{D\H) 



(32) 



where C{D\~fi ,1-1) is the likelihood (the probability of 
measuring the data, given the model "H with parame- 
ters 7^), 7r(7t|7{) is the prior (any constraints already 
existing on the model parameters) and finally p{D\'H) is 
the evidence (this is important in model selection, but in 
the subsequent analysis in this paper can be ignored as 
a normalization constant). 

Markov chain Monte Carlo (MCMC) techniques pro- 
vide an efficient way to explore the model-parameter 
space. An initial point, Xq, is drawn from the prior dis- 
tribution and then at each subsequent iteration, i, a new 
point, 1^ , is drawn from a proposal distribution, g(7/|~a^) 
and the Metropolis-Hastings ratio evaluated. 



R 



nit)C{D\-t,n)qi-^\-t) 
T:{^)C{D\-^„H)q{-^\^y 



(33) 



A random sample is drawn from a uniform distribution, 
u e C/[0, 1], and if u < i? the move to the new point is 



accepted and we set x^+i 
rejected and we set a7+i = 



. If u > i?, the move is 
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The MCMC samples can be used to carry out integrals 
over the posterior 

1 ^ 

(34) 

i—l 

The ID marginalized posterior probability distributions 
in individual model parameter follows by binning the 
chain samples in that parameter. 

The trick to using this technique efficiently is to choose 
an appropriate proposal distribution. In our analysis we 
employ an adaptive MCMC procedure, which utilises an 
"in-flight" estimation of the sampled chain's covariance 
matrix to construct an updating proposal distribution. 
This covariance matrix is updated at each iteration, with 
a certain chain memory 173^751 . We use several of the 
procedures outlined in Ref. [75|. 

For the first 100 points in the chain, simple Gaussian 
proposal distributions for each individual parameter are 
used. These first 100 points are merely used to provide 
a starting point for the covariance matrix evaluation and 
so the exact proposal distribution used in this stage is 
not important. After the first hundred points are sam- 
pled, we begin generating points via the adaptive proce- 
dure. For a £)-dimensional target posterior distribution, 
we suppose that at the i*'' iteration we have sampled at 
least H points, where the fixed integer H is the mem- 
ory parameter. We then generate a i'-dimensional vec- 
tor of trial parameters, y, via a linear mapping of an 
i?- dimensional vector of unit- variance Gaussian random 
scalars, ^ , 

t = C'/'^, (35) 

where C ' is the positive-definite square root of the 
D X D covariance matrix evaluated using the previous H 
points. The covariance matrix may be calculated by col- 
lecting the previous H points in the chain into a.n H x D 
matrix K, with each row representing one sampled point. 
Then, 

C = ;^k'^K, (36) 

where the centered matrix, K, is constructed by centering 
each column of K around the means of the respective 
parameters, calculated from the H samples. 

We then generate the trial parameter vector ~j/ via 

f ^ AA(1^, clC) ^ 1^ + -^^K^'t , (37) 

where Cd is a variable which depends only on the dimen- 
sionality of the target distribution. This variable is used 
to optimise the efficiency of the sampling process, and we 
use the value of « 2A/Vd 

With a memory parameter which is less than the total 
past history of the chain, this is denoted as the Adaptive 
Proposal (AP) algorithm [t^. Since the proposal distri- 
bution is updated constantly and relies on previous chain 



information, this procedure is not Markovian, and does 
not have the correct ergodicity properties for an MCMC 
algorithm ^73*] . In principle this can bias the reconstruc- 
tion of the target posterior; however this bias is ignorable 
in many practical applications, and for well-behaved tar- 
get posterior distributions [tsI . [t^ . If the entire previous 
chain is used to update the covariance matrix, then this 
algorithm is known as the Adaptive Metropolis (AM) al- 
gorithm [t^. The AM algorithm does not suiTer from 
the biases which can occur in the AP algorithm, and er- 
godicity is retained.'' We use the AM algorithm in our 
work. 



VI. RESULTS 

A. Posterior recovery 

An analysis using the full 10^ observations expected 
in a year of ET data is computationally prohibitive, so 
we use a working reference sample of ~ 4500 detections 
(corresponding to a shorter observation time or a lower 
merger rate density) and extrapolate to the expected 
number of detections, as discussed in Sec. I VI Cl For each 
analysis, we ran 120 independent adaptive MCMC chains 
of 5000 points on the same data catalog. We then used 
the last point from each chain to initialize a follow-up 
run of another 5000 iterations. The first 2000 points from 
each chain of the follow-up run were discarded as burn-in. 
This procedure therefore generated 360, 000 points, with 
an average acceptance rate of ~ 30%. The analysis of the 
4500-event reference catalog took ^3.5 hrs in total. Our 
sampled points were analyzed using the CosmoloGUI 
package |77| . 



B. Marginalized posterior distributions 

In Fig. 13] we show the recovered marginalized 2D pos- 
terior distributions (with 68% and 95% confidence con- 
tours) for the reference catalog. In Fig. |^c)]we observe 
a correlation between the recovered dark-energy param- 
eters. This is easily explained by the fact that a given 
cataloged luminosity distance may be consistent^ with 
a set of Wo and Wa combinations, which will depend on 
the redshift of the source. Since the majority of detected 



^ In the AM algorithm the covariance of the proposal distribution 
is actually taken to be C -I- el^ , where is the D-dimensional 
identity matrix. Choosing e > allows for the correct ergod- 
icity properties of an MCMC algorithm to be retained, and in 
practice is useful if the covariance of the chain has a tendency to 
degenerate. However, this parameter can be set very small with 
respect to the size of the target space, and in practice can be set 
to zero. 

* Here, by "consistent" we mean within ±1% of the reference value. 
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systems will be centred around z I, the wq — Wa corre- 
lation will be dominated by these sources.^ In Fig. [ ^a)| 
a negative correlation is observed between the recovered 
values for wq and /ins- For a given cataloged luminos- 
ity distance and fixed Wa, a low value of wq will imply 
a low redshift in that model. When this redshift is used 
to compute A4 from A^^, we obtain a large value of the 
chirp mass, which is consistent with a chirp-mass dis- 
tribution (and hence a NS -mass distribution) centered at 
larger values. Figure [Wb)] merely shows the combined in- 
formation of Figures [ ^a)| and [ ^c) (where the recovered 
Wa values are negatively correlated with the w p values); 
therefore the correlation observed in Fig, [ ^b)] is positive. 

A strong positive correlation is observed between the 
SFR-de nsity SF2 ansat z pa rameters, /3i and (32, as seen 
in Fig. gi^d)! while Fig. [ ^e)| shows a weak negative corre- 
lation between a and These correlations correspond 
to keeping the merger-rate density approximately con- 
stant. We calculated which combinations of a, /3i and 
were consistent with a given merger-rate density, at 
a variety of redshifts. We found that there was a strong 
positive correlation in these points between /3i and /32, 
but the correlation between a and /3i changed sign as 
the redshift increased. The greatest change occurred as 
the redshift was increased from to 1, where the corre- 
lation then reversed; however at z = 4 the magnitude of 
the correlation was still not as large as it was at z = 0.1. 
This leads us to believe that although the Dl distribu- 
tion of detected sources is peaked around ~ 6 Gpc, with 
a long tail to ^ 45 Gpc, the lower distance sources dom- 
inate the a — /?! correlation, giving an overall negative 
correlation. 

In Fig. [4l we show the marginalized ID posterior dis- 
tributions for the model parameters. The dotted lines in 
the plots indicate the 68% and 95% confidence regions of 
the marginalized distributions.^*^ 



C. Precision scaling with number of detections 

We performed similar analyses on catalogs containing 
various numbers of detections, culminating in a run with 



TABLE III: 95% confidence intervals obtained from a catalog 
of 10'^ detections, with reference parameters used to generate the 
data. AX gives the width of the 95% confidence interval. 



® The correlation between the two dark-energy EOS parameters 
can be reduced by rebinning the MCMC samples using the Wang 
parametrization [7g| . This simply involves a tranformation from 
the (wo,Wa) parametrisation to (w'Oi '"'0.5)1 where uiq.s = wq "I" 

{Wa/3). 

While these results were computed using the fast merger-rate ap- 
proximation, we also analysed a catalog using the full merger-rate 
density. The 95% confidence intervals of the marginalized pos- 
terior distributions were consistent with our approximate analy- 
sis, justifying the use of the approximation to compute the rest 
of our results. No correlations between the merger-rate density 
parameters and the dark-energy EOS parameters were found, 
which supports our earlier statement that the dependence of the 
merger-rate density on the underlying cosmological parameters 
is weak within the applied priors. 



Parameter Reference value 


95% conf. interval 


AX 




0.06 


[0.059688 , 0.060254] 


0.000566 




1.35 


[1.347408 , 1.351789] 


0.00438 


Wo 


-1.0 


[-1.036403 , -0.949623] 


0.0869 


Wa 


0.0 


[-0.195630 , 0.073602] 


0.269 


a 


-1.0 


[-1.026691 , -0.961659] 


0.0650 


/3i 


3.4 


[3.318136 , 3.605810] 


0.288 


/?2 


3.4 


[3.310287 , 3.582895] 


0.273 



10^ detections. We can characterize the precision with 
which we can measure the various model parameters by 
the 95% confidence intervals. Recording these intervals 
for all parameters for varying catalog sizes, and dividing 
by the reference sample intervals gave the results shown 
in Fig. [S] This clearly shows that the precisions scale 
as 1/ViVo as we would expect. Parameter measurement 
accuracies for the 10^-event catalog are shown in Table 
mil We see that the measurement precisions of the dark- 
energy EOS parameters are the same order of magnitude 
as those forecast for CMB+BAO+SNIa [13, as discussed 
in Sec. IVD] 



D. Including and accounting for errors 

Distance measurements from a third-generation GW- 
interferometer network will not be error-free. While a 
network consisting of a single FT plus one other right- 
angle interferometer can place constraints on a source's 
sky location and luminosity distance, the precisions of 
these properties are improved to almost the 3-ET net- 
work level by the inclusion of a second additional right- 
angle interferometer [sHj- The redshifted chirp mass is 
expected to be very well constrained (< 0.5% error J4(f]), 
and so we ignore measurement errors in this parameter. 
We assume the error in the luminosity distance arising 
from instrumental noise scales as ~ 1/p, and include the 
effects of weak lensing as a further source of error. The 
weak-lensing error on luminosity distance measurements 
at z ~ 1 is approximately 5%, and we linearly extrapo- 
late this to all other redshifts jll, [H [zi, IsS] • While sev- 
eral techniques have been proposed to reduce this weak- 
lensing error [e.g. Refs. [Sll i82] and references therein], 
we assume no correction has been done, corresponding 
to a worst-case scenario. 

Errors on the distance-redshift relation from binary- 
system peculiar velocities are much smaller than instru- 
mental and weak-lensing errors at all but the lowest red- 
shifts, becoming comparable with these at z ~ 0.1 where 
the error is < 1%, and decreasing sharply at higher red- 
shifts [Ref. i83j] and references therein]. The lowest red- 
shift in our reference catalog is ~ 0.05, where the pe- 
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FIG. 3: Marginalized 2D posterior distributions for the reference catalog of 4500 detections. Only those 2D distributions showing 
correlations between parameters are shown. The reference parameters are fj,^Q = 1.35Mq, CTfjS = 0.06Mq, wq = — 1, Wa = 0, a = — 1 
and /3i = /32 = 3.4. 




FIG. 4: Marginalized ID posterior distributions for the reference catalog of 4500 detections. Dotted lines indicate the boundaries of 
the 95% and 68% confidence intervals. The reference parameters are /^ns = 1.35M0, cr^g = O.OOMq, wq = —1, Wa = 0, a = —1 and 
/3i = 3.4. 



culiar velocity errors will dominate, but only lead to an 
error of < 2%. The sensitivity of the luminosity distance 
to the dark-energy EOS parameters is very weak in this 
redshift regime; hence peculiar velocities are unlikely to 
introduce significant parameter bias/inaccuracy, and we 
ignore them here. 

We also ignore the effect of detector calibration er- 
rors, which, unlike statistical measurement uncertainties, 
would not be mitigated by boosting the detection rate. 
Such systematic biases have recently been studied in the 
case of advanced-era detectors [111 , and found to induce 



a systematic shift in the estimated system parameters 
which is a small fraction of the statistical measurement 
errors. We ignore waveform-modelling errors in our anal- 
ysis, since current post-Newtonian models will only break 
down close to the onset of the merger-phase, and for the 
neutron-star binaries considered in this analysis this is 
at frequencies where the instrumental noise is high and 
which therefore do not contribute much to the overall 
signal-to-noise of the system. Furthermore, luminosity- 
distance determination comes primarily from the network 
triangulation which will not be significantly affected by 
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TABLE IV: 95% confidence intervals derived from the reference sample (4500 detections), both in the case where distances are 
measured precisely, and when distance errors are included and accounted for (using the error averaging technique described in the text, 
with various numbers of points sampled from the distance posterior PDF). AX gives the width of the 95% confidence interval. 



Parameter 


No errors 


Errors (50 pts) 


Errors (100 pts) 


Errors (400 pt 


^) 




95% conf 


interval 


95% conf 


interval 




95% conf 


interval 




95% eonf 


interval 


AX/AX„f 


o-Ns/M© 


[0.058785 , 


0.061447] 


[0.066378 , 


0.071911] 


2.07851 


[0.063815 , 


0.069409] 


2.10143 


[0.059309 , 


0.064849] 


2.08114 


MNs/M0 


[1.339198 , 


1.358745] 


[1.329060 , 


1.354066] 


1.27928 


[1.335499 , 


1.361690] 


1.3399 


[1.335782 , 


1.359339] 


1.20515 


Wq 


[-1.145894 , 


-0.749671] 


[-0.880092 , 


-0.338642] 


1.36653 


[-1.146052 , 
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FIG. 5: 95% confidence intervals of the ID marginalized 
distributions relative to those of the 4500-event reference catalog, 
shown as a function of the number of cataloged events. The same 
intrinsic parameters of the underlying distributions are used to 
create the mock catalogs. The expected ~ l/y/N relationship is 
overlaid on the plot. 



modelling uncertainties, and so the distance error will 
be dominated by instrumental- noise and weak-lensing, 
as discussed earlier. Similarly, the error in the distri- 
bution of possible source-redshifts arising from the mea- 
sured redshifted chirp mass will be dominated by the 
intrinsic width of the NS mass distribution rather than 
the small error in the redshifted chirp mass coming from 
instrumental noise and modelling uncertainties. 

We repeat the run of the working reference sample, off- 
setting the cataloged luminosity distance by an amount 
drawn from a Gaussian distribution, with mean at the 
true distance, and standard deviation. 



(38) 



The data collected for each event will actually be in the 
form of posterior probability density functions (PDFs) 
for the parameters, where previously we have assumed 
these are 5-functions at the true values. If the offset 
luminosity distances are assumed to be the true values 
with a 5-function posterior PDF, the chain does not move 



away from it's starting point. Hence we must take these 
errors into account in the likelihood calculation stage. 

We can account for these errors in the analysis, by 
modifying the likelihood in Eq. ([29]) to 



p I it = ~^ 



'n^(Ati7t) 



i=l 



d'^'Aid'^Aa-.-d'^A^, (39) 



where, in our case, the number of parameters k — 2, and 
~^ is the detector output, which is a combination of No 
signals, /li, and noise, it . Equation pop is an integral 
over all possible values of the source parameters that are 
consistent with the data. The first term inside the square 
bracket is the computed posterior PDF for the detected 
population of sources. Concern has been raised that the 
high event rate of ET detections may lead to a confusion 
background. However, the noise-power- spectral- density- 
weighted signals are short enough that there is not ex- 
pected to be significant overlap [161 |4Q|. Hence, these 
detections should be uncorrelated, with independent pa- 
rameter estimates, and so this first term reduces to the 
product of the posterior PDFs for each detection. 

If the posterior PDF for a given source has been 
obtained via IMCMC techniques, then the integral in 
Eq. (|39p may be computed by summing over the chain 
samples. Thus, errors may be accounted for by making 
the following replacement in Eq. ([29]) : 



* .7 = 1 



(40) 



where Mi is the number of points in the chain for the i*'' 
source's PDF, and A-'^^ is the element of the discrete 
chain representing this PDF. This technique does not 
assume a specific form for the PDF, and can be used 
in the case of multimodal distributions. 

We represent the Dl posterior PDF for each source 
by a chain of 50 points, drawn from a normal distribu- 
tion with standard deviation as in Eq. (j38p . and a mean 
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equal to the value in the data catalog, which in this anal- 
ysis, as discussed earlier, includes an error to offset it 
from the true value. Results are shown in Table IIVI We 
see that a significant bias in the reconstructed model pa- 
rameters still exists. We suspected that this bias arose 
from using only 50 points to evaluate the distance pos- 
terior PDFs. We therefore repeated the analysis with 
an increasing number of points sampled from the dis- 
tance posterior PDF.^^ With 100 points, all biases are 
corrected expect for that in cnS: and the ratio of the 
95% confidence interval widths to the reference widths is 
not significantly different from the 50 point case. This 
suggests that a larger number of points in the error av- 
eraging technique will be necessary to correct all biases, 
but this is not neccesary to estimate parameter measure- 
ment accuracies in the presence of distance errors. Using 
400 points sampled from the distance posterior PDF all 
biases in the parameter posterior distributions appear to 
be corrected, in the sense that the reference parameters 
then lay within the 95% confidence intervals of the ID 
marginalized posterior distributions. 

Overall, we find that the result of properly account- 
ing for instrumental and weak-lensing errors is that pa- 
rameter measurement precisions are, at worst, approxi- 
mately halved. Thus instrumental and weak-lensing in- 
duced errors should not affect our general conclusions 
about the science capabilities of a third-generation GW- 
interferometer network. We carry out the remainder of 
this study using catalogs which are generated and ana- 
lyzed without including errors. 



E. Precision scaling with intrinsic parameters 

We now investigate how the ability of ET to constrain 
the parameters of the underlying distributions is affected 
by the values of the intrinsic parameters themselves. This 
is similar to the kind of analysis performed in our pre- 
vious study with second-generation interferometers [22l |. 
We perform analyses of data catalogs generated with dif- 
ferent intrinsic parameter combinations; multiple runs 
are performed on each parameter combination. We vary 
one parameter at a time, fixing all others to the reference 
values. 

Varying the intrinsic parameters with fixed SNR 
threshold will alter the expected detection rate. This 
is illustrated in Table |Vl where the model with reference 
parameters has an expected detection rate of ~ 10^ yr~^. 
The expected ^ 1/\/N relationship is well established, 
as shown in Fig. [5] Hence, we remove this number effect 



The posterior distributions obtained via this analysis should be 
considered estimates of the true distributions, since the long like- 
lihood computation time required by this error-analysis means 
that we did not collect as many samples as when errors are ig- 
nored. We performed burn-in runs, and then follow-up runs to 
estimate the posterior distributions as well as feasible. 



by generating catalogs with the same number of events 
(4500 each in order to compare against the reference cat- 
alog) . Therefore we are testing how the cosmological, as- 
trophysical and intrinsic-mass distributions of coalescing 
DNS binaries impact the precision of parameter recovery. 

The results of these analyses are shown in Fig. [S) We 
see that as uns is increased the measurement precision of 
both the NS mass distribution and dark-energy EOS pa- 
rameters decreases. We found a similar trend in Ref. [23]. 
This makes sense, since if we have an intrinsically narrow 
NS mass distribution, then we have a good idea of what 
the intrinsic masses of the systems are and the range of 
candidate redshifts produced from a measured redshifted 
chirp mass will be narrow, improving the precision with 
which we can recover cosmological parameters. 

A variation in the intrinsic /iNS (not shown) pro- 
duces accuracies comparable to the reference accuracies. 
Hence, the impact of the intrinsic value of the NS mass- 
distribution mean on parameter accuracies is predomi- 
nantly through the change to the expected detection rates 
i.e., a larger mean, on average, will lead to larger chirp 
masses, so that detections can be made from a larger 
volume (see Eq. ([5])). 

Increasing the value of the EOS parameter wq in- 
creases the precision with which we can recover wq, Wa 
and /iNS- As wq is increased, while the intrinsic Wa is 
fixed at zero, the recovered posteriors for these param- 
eters are squeezed by the prior restrictions, wq < —1/3 
and {wq + Wa) < —1/3. A larger intrinsic wq increases 
the horizon distance of detections, which permits greater 
sensitivity to the dark-energy EOS parameters. Fur- 
thermore, a narrowed range of cosmological parameters 
means that the range of candidate redshifts is also nar- 
rowed, such that the precision of the recovered NS mass- 
distribution mean (deduced from the redshifted chirp 
mass) improves. We also notice these effects when the 
intrinsic Wa is increased with the intrinsic wq fixed at the 
reference value. However, the effect is less pronounced in 
this case, since Wa is a first-order correction to Wq. 

As the power-law index, a, is increased the average 
delay between the formation of the massive progenitor 
system and the merging of the final compact-system in- 
creases. This means that more systems formed at higher 
redshifts survive to merge at lower redshifts, and hence 
the merger-rate density is boosted to higher values at 
lower redshifts. In addition, as a increases the merger- 
rate density tracks the underlying SFR-density to a lesser 
extent, so it becomes more difficult to extract the details 
of the SFR-density. Hence the /3i.2 distributions widen 
to reflect this reduced sensitivity to the underlying SFR- 
density. 

When we keep the intrinsic values of /3i and /32 equal 
(not shown), we find that varying these by ±0.4 around 
the reference value has a negligible impact on the mea- 
surement precision of the parameters. A higher common 
/3i.2 value leads to a larger expected detection rate, but 
this is a small effect. 

Lowering the intrinsic value of with /32 fixed, shifts 
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TABLE V: The expected detection rates for different choices of intrinsic parameters of the underlying distribution are compared to 
the reference expected detection rate. One parameter is varied at a time, with all other parameters kept fixed at their reference values. 
The variation of the expected detection rate with the intrinsic value of ctns is not shown, since this parameter is not used in the rate 
calculation (see Appendix [Bll . 
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FIG. 6: The variation of measurement precision with different choices of the intrinsic parameters of the underlying distributions. One 
parameter is varied at a time, and in the interest of testing how the precision of parameter recovery is affected by the underlying 
distribution of events, all catalogs are generated with the same number of events (4500 to match the reference catalog). Each point in 
each panel represents the average 95% confidence interval width of 3 realizations of the catalog. 



the distribution of events to lower distances, and changes 
the shape of the underlying merger-rate density. This 
distribution is consistent with a wider range of a vahies 
than the reference distribution, since the sensitivity of 
the merger-rate density to a is reduced at lower redshifts. 
This causes the marginahzed a-posterior distribution to 
widen. The same is true when the intrinsic value of is 
increased. 



F. Varying the SNR threshold 

We also generated catalogs for different SNR thresh- 
olds defining the detectability of merging systems. Mul- 
tiple catalogs were analyzed for each SNR threshold, but 



once again the number of events was fixed at 4500 to 
match the reference catalog (see Fig. [71). An increase in 
the SNR threshold is equivalent to the characteristic dis- 
tance reach of the detectors decreasing. Hence we would 
expect the sensitivity of the data to varying dark-energy 
EOS parameters, which have a greater influence at larger 
redshifts, to be reduced. A reduced characteristic reach 
would also result from a larger low- frequency cutoff, in 
the detector's noise power spectral density. In the recent 
mock ET data challenge, Regimbau et al. [i^l found that 
confusion between two or more signals rarely affected the 
analysis performance when /; = 25 Hz. Standard algo- 
rithms currently employed for LSC- Virgo analyses can- 
not handle templates longer than a few minutes; however 
multiband filter methods are being developed which will 
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FIG. 7: We show the variation of the expected detection rate as 
the SNR threshold, po, is raised. This can also be interpreted as 
lowering the characteristic distance reach of the network. Since p 
scales as 1/Dl (see Eq. 10), and the difference between the 
luminosity distance and radial comoving distance becomes smaller 
at lower redshift, one would expect that at high enough values of 
po the comoving detection volume (and hence detection rate) 
would scale as l/pg. This is approximately valid for po > 15. 



allow /; to be pushed below 25 Hz. In Fig. |8] one can 
see that with // = 25 Hz, the effective SNR threshold is 
raised from the reference value of 8 (with fi — 1 Hz) to 
~ 12.4. 

From Fig. [ ^a)[ we see that as the SNR threshold is 
increased, with the number of cataloged events fixed, the 
accuracies of /3i and (32 degrade sharply. At higher SNR 
thresholds (or, equivalently, at lower distance reaches) 
the sensitivity of the merger-rate density to varying /3i_2 
is reduced; hence the wider posterior distributions. The 
measurement precision of a increases slightly as the SNR 
threshold is increased from 6 — 12. One might expect a 
to show the same trend as /3i and /32 , since an increasing 
SNR threshold pushes the events to lower redshifts where 
the sensitivity of the merger-rate density to a is reduced. 
However, we see in Fig. [10] that the merger-rate density, 
for various choices of a (but all other parameters fixed), 
is relatively featureless beyond ~ 2.5. The distribution 
of the merger-rate density in the redshift window of 
2.5 — 7 could be approximately linearly scaled to satisfy 
a large range of a}'^ Therefore, given that our likelihood 
statistic is insensitive to linear scalings of the merger-rate 
density [see Eq. the significant number of high- 

redshift detections in a po = 6 catalog will widen the a 
posterior distribution, while most a-information is found 
in the redshift window ~ 1 — 2, where the merger-rate 
density has more features. 



The same argument did not apply when /3i^2 was varied, since 
this not only shifted the distribution to lower redshifts but al- 
tered the shape of the merger-rate density in a way that could 
not be equated with a linear scaling in any redshift window. 



In Fig. [ ^b)[ we see that the measurement accuracy 
of wq and Wa is slightly reduced for higher SNR thresh- 
olds; this is a small effect, and is expected with a cat- 
alog shifted to lower redshifts, where distances are less 
sensitive to varying dark-energy EOS parameters. The 
accuracy of wq only varies by ~ ±5%, since we remain 
sensitive to detections at tens of Gpcs even with an SNR 
threshold of 12. However, Wa shows a stronger variation 
since it is a higher order correction to the EOS parame- 
ter, and distances become sensitive to this parameter at 
higher redshifts than they do to wq. We also see that 
the measurement precision of /xns and ctns is increased 
slightly as we move to larger SNR thresholds. This small 
effect is probably due to the fact that a lower redshift 
range in the data catalog will mean that the redshifted 
chirp mass is closer to the intrinsic chirp mass. 

Although this suggests that a greater distance reach 
will improve the precision of cosmological parameter re- 
covery, we have so far ignored distance errors. In fact, 
instrumental and weak-lensing errors impart an interest- 
ing redshift evolution to the WQ-sensitivity, which we ap- 
proximate as jSGj l 



Aw{z) 



dwo 



X DlX v/(l/p)2 + (0.05z)2. (41) 



In Fig. [TT] we see that the sensitivity of the luminosity 
distance to the cosmological parameter wq is greatest at 
z ~ 1, since wq has a very weak intrinsic impact on Z?^ at 
low redshifts and distance errors dominate at higher red- 
shifts. Increasing a detector's distance reach will raise the 
fraction of high-redshift cataloged events. We calculate 
the effective measurement precision of wq from our refer- 
ence catalog by adding the Awq values from each event 
in quadrature i.e. l/Ai(Jo,cff = -\/^(l/Ait;o,i)^- This 
is repeated for various lower and higher SNR threshold 
values. We perform these calculations for catalogs con- 
taining the same number of events (4500 to match the 
reference catalog), and for catalogs with the number of 
events scaled by the ratio of the expected detection rate 
for each SNR threshold to the reference threshold (which 
in this analysis is 8). The results are shown in Table 
IVIl where we see that for catalogs with the same num- 
ber of events, lowering the SNR threshold actually wors- 
ens the precision of wq recovery since the distribution of 
events is weighted to higher redshifts, where distance er- 
rors degrade the precision. Increasing the SNR threshold 
reduces the number of events at high redshift and hence 
mitigates the degradation of precision due to distance er- 
rors (see Fig. [TT|) . However, this effect slows down with 
increasing SNR threshold. For catalogs with numbers of 
events scaled to match the expected detection rate for 
each SNR threshold, we see that the increased number of 
events associated with a lower SNR threshold is enough 
to compensate for degradation of precision from higher 
redshift events. However this loss of precision means that 
lowering the SNR threshold does not lead to the 1/^/N, 

or 1//9q^^ improvement in parameter measurement preci- 
sion which one would naively expect. 
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FIG. 8: The reduction of the characteristic distance reach associated with raising the low-frequency cutoff, of 3'^'^-generation 
detectors. This can also be interpreted, via Eq. I I28I I. as raising the network SNR threshold. The figures were produced using the ET-D 
noise curve [35|| . 
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FIG. 9: The variation of the parameter measurement precisions with the network SNR threshold. The left panel shows precisions, 
characterized by the width of the 68% confidence intervals, for the merger-rate density parameters, while the right panel shows 
precisions, characterized by the 95% confidence intervals used elsewhere, for all other parameters. We use the narrower confidence 
intervals for the merger-rate parameters to mitigate the effct of poor sampling in the low-a region which was observed in some 
AM-MCMC chains in this analysis. All catalogs contain the same numbers of events at each threshold value, which, as in the previous 
subsection, is 4500 to match the reference catalog. 



TABLE VI: The events from catalogs with different SNR 
thresholds are used to compute an effective uiq precision, by 
adding the Auiq values of each event in quadrature. This analysis 
is performed for catalogs with the same number of events, and for 
catalogs with the numbers of events scaled to match the expected 
detection rate for each SNR threshold. 



A-mo,cff7lO 
(jV„ = Af,ef) (JVq = / X Af,,f) 



12 
20 
30 



1.64 
1.00 
0.399 
0.0936 
0.0271 



8.08 
6.82 
5.33 
4.00 
3.38 



6.33 
6.82 
8.35 
13.1 
20.6 



Finally we address the issue of having assumed that 
Earth motion does not modulate the antenna patterns 
of the detectors. The time spent "in-band" by an 
inspiraling-event scales as [1^ 



5.4 



1.22Mp 



-5/3 



f,. 



-8/3 



days. 



(42) 



Hence, a detector with a low-frequency cutoff of 1 Hz 
(as we have assumed) could have events in band for as 
long as ~ 5 days. In this case, a correct treatment of the 
antenna pattern modulation would be needed. However, 
if we increase /; to ~ 8 Hz, then the maximum time 
spent in-band is less than 30 minutes, and ignoring the 
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FIG. 10: We show the redshift distribution of the DNS 
merger-rate density for various choices of the power-law index of 
the delay-time distribution, a (all other parameters are fixed at 
their reference values). The merger-rate density is relatively 
featureless beyond z ~ 2.5, making it difficult for our analysis 
(which is insensitive to linear scalings of the merger-rate density) 
to discriminate between values of a in this range. Overlaid on 
this figure we show the redshift distribution of detections for 
various SNR thresholds. 
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FIG. 11: The redshift dependence of the sensitivity of the 
luminosity distance to wg . This parameter has a very weak 
intrinsic impact on at low redshifts, whilst distance errors 
from instrumental noise and weak-lensing dominate at higher 
redshifts. This results in a redshift "sweet-spot", where these 
effects are minimised for lines of constant SNR. We also plot the 
individual Awq values calculated for the reference catalog events. 



antenna pattern modulation is reasonable. In Fig. [8] we 
see that a low-frequency cutoff of 8 Hz is equivalent to 
raising the SNR threshold to ~ 9, and from Fig. ^ we 
see that the precision of parameter recovery is within ~ 
±10% of the reference precisions for an SNR threshold of 
9. Therefore our approximate treatment of the network 
antenna patterns would seem reasonable. 



VII. CONCLUSIONS 

We have built on our previous work j22| which explored 
the capabilities of an advanced (i.e. second-generation) 
GW-interferometer network to constrain aspects of the 
NS mass distribution in DNS systems, as well as the 
Hubble constant. The technique we employed used only 
information obtained via analysis of the GWs detected 
in such a network. In this paper we extended the anal- 
ysis to a possible third-generation network, consisting 
of the proposed Einstein Telescope, and complemented 
by third-generation right-angled interferometers at LIGO 
Livingston and LIGO-India. The target for the Einstein 
Telescope is a broadband factor of 10 sensitivity increase 
with respect to advanced detectors, but to also extend 
the low-frequency sensitivity of ground-based GW inter- 
ferometers below 10 Hz. The current design for a sin- 
gle ET consists of three overlapping interferometers, ar- 
ranged in an equilateral configuration with arm-opening 
angles of 60° [26l - l28| . Each interferometer will consist of 
a cryogenically-cooled, underground low-frequency detec- 
tor, and a high laser-power, high-frequency detector in a 
"xylophone" configuration - these two detectors work in 
tandem to suppress noise over the entire band 29, 34]. 
Current projections for funding and construction of ET 
place "first-light" sometime in the mid-2020's. 

The sources of interest in this paper are inspiraling 
double NS systems, which could be observed at rates 
of ~ 40 yr~^ by advanced detectors [ll| and rates of 
0(10^ — lOf) yr~^ may be achieved by a third-generation 
network [la, HI, E^l • These sources are commonly re- 
ferred to as self-calibrating standard sirens, since their 
distance from us is directly encoded in the emitted GWs. 
Combined with a method of redshift determination, these 
sources can be used to probe the distance-redshift rela- 
tion and hence extract constraints on background cosmo- 
logical parameters which are independent of the cosmic 
distance ladder 

Our method of cosmography using only GWs relies 
on the narrowness of the distributions of masses of NSs 
in these DNS systems. Recent analysis indicates that 
this mass distribution is indeed narrow, with a Gaussian 
mean of ~ 1.35M0 and standard deviation of O.O6M0, 
which i nay be a product of a distinct evolutionary path 
[lEIUjllS- Using a measurement of a source's redshifted 
chirp mass, we can therefore obtain a narrow candidate 
redshift distribution. A narrower intrinsic NS mass dis- 
tribution will obviously mean the precision of redshift de- 
termination increases. We can combine these with GW- 
interferometer network determinations of the luminosity 
distance to constrain cosmological parameters. 

We used a Bayesian theoretical framework to assess the 
capability of a third-generation network to measure cos- 
mological and astrophysical parameters. We performed 
7-dimensional adaptive MCMC analysis on the catalogs 
of detections, using reference parameters Hq — 70.4 
kms-iMpc-\ flmfi = 0.2726, Qkfi = -0.0006, wq = -1, 
Wa — 0, /^Ns = I.35M0 and a^s = O.O6Af0. Keeping Hq, 
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flm.o and f^A^o fixed, we found that the measurement pre- 
cisions of the dark-energy EOS parameters possible with 
a 10^-event catalog were of the same order of magnitude 
as forecasted constraints from future CMB+BAO+SNIa 
measurements [l3|- Furthermore the power-law index of 
the merger delay-time distribution, a, and the parame- 
ters of the underlying star-formation-ratc (SFR) density 
were constrained to within ~ 10%. Accounting for mea- 
surement errors degraded precisions by a factor of < 2, 
while increasing the network SNR threshold required for 
detection from 8 to 9 (which is equivalent to consider- 
ing only the ~ 30 minute section of inspiral above 8 Hz) 
changed the precisions by only ^ 10%. 

We also investigated how the precision of parameter 
recovery scaled with the values of the intrinsic parame- 
ters themselves, keeping the number of detected events 
fixed to factor out pure number-of-event effects. Varying 
the intrinsic ctns showed a linear scaling of parameter 
precision, with narrower intrinsic NS mass distributions 
favouring tighter parameter constraints. The precisions 
of the merger-rate density parameters did not appear to 
be affected in this case. Increasing the intrinsic wq and 
Wa had the effect of increasing their measurement pre- 
cision, as well as that of ^ns- This was probably due 
to the fact that larger wq and Wa give detections out to 
greater distances, where the sensitivity to these param- 
eters is higher. Tighter cosmological constraints implies 
narrower candidate redshift distributions from the cata- 
loged distances, which improves /xns precision. Increas- 
ing the intrinsic value of a meant that the merger-rate 
density tracked the underlying SFR density to a lesser 
extent and hence worsened the precision of SFR-density 
parameter recovery. As we changed the shape of the un- 
derlying SFR density to favor closer detections, the mea- 
surement precision of a worsened, since the sensitivity of 
the merger-rate density to a is lower at lower redshifts. 

Finally, we varied the criterion for a network detection, 
which we denoted by a threshold value of the network 
signal-to-noise ratio. This could also be interpreted as 
varying the characteristic distance reach of the network, 
which, in turn, could be caused by varying the detec- 
tor's low-frequency cutoff. Varying the SNR threshold 
between 6—12 caused a slight decrease in wq and Wa 
precision, as catalogs with lower distance events will be 
less sensitive to these cosmological parameters. However, 
catalogs with, on average, closer events will provide bet- 
ter NS mass-distribution parameter precision, since the 
redshifted chirp mass will be less offset with respect to 
the intrinsic chirp mass. Increasing the SNR threshold, 
and hence decreasing the characteristic distance reach of 
the network, caused a significant decrease in SFR-density 
parameter precision, since the merger-rate density is less 
sensitive to the SFR-density parameters at lower redshift. 

While the sensitivity of distances to the dark-energy 
EOS will obviously be intrinsically weak at low redshifts, 
distance-measurement errors begin to dominate at higher 
redshifts. We found that for a fixed number of events in 
a catalog, lowering the SNR threshold actually worsened 



the precision of wq recovery since events are weighted to 
higher redshifts, where distance errors degrade the mea- 
surement precision. The larger expected detection rate 
associated with lower SNR thresholds is enough to re- 
verse this effect, but means that lowering the SNR thresh- 
old (or increasing the network's distance reach) does not 
lead to the great improvement in parameter measurement 
precision which one would naively expect. 

We have not considered association of GW detections 
with an EM counterpart, either through precise sGRB 
[H, [13 or host-galaxy association. The latter technique 
may only be possible with ~ 0.01% of detectable GW 
events ^J]. However, the technique we have used has 
been shown previously [12] to be well complemented by 
precision redshift information. In particular, we found 
that if redshift information (measured to much greater 
precision than the luminosity distance) is available for ^ 
10% of a GW-event catalog, then measurement precisions 
of parameters were more than doubled. 

This paper completes our proof-of-principle study of 
this GW-only cosmographic technique. We have shown 
the significant potential for a third-generation network 
including the Einstein Telescope to place interesting con- 
straints on the NS mass distributions in DNS systems, 
the dark-energy EOS, the average delay between the 
formation of the DNS-system progenitors and the final 
merger and the underlying SFR density in the Universe. 
Over the following decade tighter constraints will be de- 
rived for the NS mass distribution, delay-time distribu- 
tion of DNS systems, and the SFR density, which can be 
readily incorporated within this technique. We intend to 
test this technique in the upcoming FT mock data chal- 
lenge, as well as study the ability of this technique to 
discriminate between NS mass distributions from differ- 
ent metallicity progenitors, different delay-time distribu- 
tions resulting from different formation paths, and possi- 
bly multimodal NS mass distributions. Unshackling GW 
cosmography from its reliance on EM counterparts will be 
an important step in establishing DNS systems as phys- 
ical distance indicators, and contribute to GW-analysis 
becoming a precision astrophysical tool. 
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Appendix A: DNS Merger-rate density 

Here, we provide a more detailed discussion of the as- 
trophysics of DNS mergers and justification for the ansatz 
we employ for the merger-rate density. 



1. Merger-delay distribution, dPm/dt 

Assuming the number, TV, of DNS binaries born with 
separation a follows dN/da cx (s^l, the merger-delay 
distribution is 



dPm 

dt 



dN 
dT„ 



dN da 
da dTur 



(Al) 



An early suggestion by Piran was to consider newly 
formed DNS binaries as having the same orbital sepa- 
ration distribution as normal-abundance main-sequence 
stars [13, [5^ . For normal-abundance main-sequence bi- 
nary systems, the distribution of periods has been found 
to be flat in In (P) (where P is the binary period) (55| . 
or lognormal [83, |8g| . If we follow Ref . 54] , and ignor- 
ing the progenitor ellipticity distribution, then the ini- 
tial DNS orbital separation distribution will be flat in 
In (a) (where a is the semi-major axis of the binary or- 
bit). Therefore 7 = -1 and dP^/di cx t'^^^ 

The catalog of DNS systems in Ref. [56| is used by 
Refs. [13, to estimate the merger-delay distribution 
of observed systems; it approximately follows (1/i), but 
there appears to be an excess of systems below an inspiral 
time of 100 Myr. Selection effects having to do with the 
difficulty of detecting binary pulsars in close orbit, due to 
the large and rapidly varying Doppler shift, may signifi- 
cantly affect the reconstructed merger-delay distribution 
below a few hundred Myr. The authors comment that 
with such a small sample of systems 6) it is difficult to 
make predictive conclusions about this distribution, but 
that it is the best one can presently do with the observa- 
tions. 

Population synthesis calculations in Refs. [13, Ull ap- 
pear to show the cumulative merger-delay distribution 
being approximately linear in In (t) (in which case the 
PDF varies as ~ t~^), while the studies in Refs. ^59n61j 
show that a (1/t) PDF is an appropriate approximation 
over several orders of magnitude of the delay-time. Fur- 
thermore, the population synthesis calculations of Ref. 



The caveats here are that eUipticity can have a significant efl'ect 
on inspiral timescales, and it is not obvious that DNS systems 
should have the same orbital separation distribution as main- 
sequence binaries, since the two supernovae the systems survive 
would likely modify it. Furthermore the distribution functions 
for progenitor evolutionary timescales and merger timescales are 
not independent. The evolutionary timescale depends on the 
mass of the progenitor system components and the gravitational 
inspiral timescale depends on the chirp mass of the system. 
Strictly speaking the joint probability consideration should be 
considered [s^; however we ignore this subtlety here. 



[62| . in their study of the formation rates of short- and 
long-GRBs, indicates an approximate (l/t) delay-time 
distribution for NS-NS and NS-black-hole systems. 

Population synthesis calculations have also proposed 
previously unconsidered DNS-formation channels; specif- 
ically Ref. [ssl suggests a formation channel via a double 
common-envelope phase between two low-mass helium 
stars, such that the subsequently formed NSs would not 
have had time to accrete matter and be recycled. These 
DNS systems would be under-represented in Galactic 
pulsar surveys, since they would be observable as radio 
pulsars for a much shorter time scale than recycled pul- 
sar systems. Hence DNS coalescence rate calculations 
based on the observed Galactic pulsar sample need to 
take into account any observational biases. Another of 
these new formation channels involves a stage of hyper- 
critical common-envelope accretion from a low-mass he- 
lium giant to the firstborn NS, resulting in a popula- 
tion of tight, short-lived binaries (with merger-timescales 
< 1.0 Myr) which may contribute significantly to the to- 
tal number of coalescences [ooj . 

In the 2004 study of merging DNS systems as the 
source of sGRBs by Ando [91|], the merger-delay distri- 
bution is modelled as a power- law (^ t"), and the cal- 
culated GRB rate densities are found to be relatively 
insensitive to the lower-cutoff time necessitated by such 
a parametrization, but considerably sensitive to a. The 
characteristic upper inspiral timescale is also of interest; 
several known DNS systems are calculated to have inspi- 
ral times exceeding 10 Gyr (Refs. [5^ H^l and references 
therein). If these are representatives of a class of DNS 
systems resulting from a different evolutionary path than 
the lower time scale systems, then this evolutionary path 
need not be considered for detectable GW sources. 

With the above considerations in mind, in the present 
study we adopt a power-law merger-delay distribution for 
DNS systems, with a reference index of —1. This is sup- 
ported by the existing (albeit sparse) observational data 
on the gravitational inspiral times of Galactic DNS sys- 
tems, and the prevalence of power-law delay distributions 
found in population synthesis studies. For normalization 
purposes, we adopt a lower delay time of 50 Myr, since 
the massive progenitor system (containing components 
with masses between ^ 8 — 20Mo for NS-NS system for- 
mation) may require an evolutionary timescale of > 50 
Myr.^"' Taking this as a lower delay time avoids con- 
siderations of the (possibly significant) DNS-formation 
channel with the extra mass-transfer episode (which cre- 
ates a peak in the delay-time distribution around ^ 20 
Myr, and corresponds to a population of tight, short- 
lived DNS systems [90]). We model only DNS systems 
formed via the classical formation channel [59, 0, |94| 



This evolutionary time scale is an approximate main-sequence 
lifetime for a 8Mq star, burning ~ 10% of its core hydro- 
gen, and obtained via the simple scaling relationship, TovoI ^ 
lO4(M/M0)-2-5 Myr. 
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for which the ^ t^^ delay-time distribution is an ap- 
propriate approximation over several orders of magni- 
tude. The power-law index will have a greater impact 
on merger-rate density calculations than the lower cutoff 
time. We assume an upper inspiral timescale equal to 
the cosmology-dependent age of the Universe. For the 
present study, the power-law index in this merger-delay 
distribution is labeled a. 



2. Star- formation rate density, dp^/dt 

The determination of the low-redshift SFR density has 
been achieved via a wide variety of techniques, utiliz- 
ing light at different wavelengths. However, these mea- 
surements become more difhcult at higher redshifts since 
many of the techniques successfully employed in the low-z 
Universe rely on light at wavelengths that cannot be de- 
tected beyond z ~ 4. This leaves us with only a handful 
of available techniques to probe the high-redshift star- 
formation history. 

The estimations of star-formation rates at high- 
redshift are obtained from measurements of UV lumi- 
nosity functions (LPs) , which tells us how many galaxies 
emit light in the UV-band in a given epoch. Excepting 
galaxies with the largest SFRs (which likely suffer from 
significant dust extinction) UV light has been shown to 
be a good tracer of the SFR (Refs. [13, IHI and refer- 
ences therein) . Dust-extinction of UV light can be inves- 
tigated, and hence corrected for, via the measurement of 
the UV-continuum slope, which has been shown to be 
well-correlated with dust extinction in the local Universe 
(Ref. i64| and references therein). A systematic study 
of the high-redshift SFR density was undertaken in Ref. 
[G^ I using Hubble Space Telescope data. Correcting their 
UV luminosity density calculations for dust extinction, 
and converting this to an estimate of the SFR density, 
yielded significant evolution of the SFR density between 
< ^ ^ 6. The SFR density is shown to rise out to 
z ~ 2 — 4, followed by a decrease out to z ~ 6. This 
decrease is shown to continue out to z ~ 8.5 [65 1. 



Given that only a handful of techniques exist to probe 
the high-redshift star-formation history, we will have to 
wait until further studies are carried out, or new tech- 
niques are developed, to complement the analyses in Refs. 
[GJ, 111] . In our present study, we are only interested in 
a sensible model of the redshift evolution of the SFR 
density, which we can parametrize for a Bayesian infer- 
ence analysis. Several of the studies (Refs. [13, [H, l9l|) 
mentioned in Sec. IIIIB 2[ as well as several other stud- 
ies which attempt to fit GRB densities to delayed SFR- 
density models (e.g., Ref. i96i]), employed the SF2 model 
of Porciani and Madau 1631 . Of the three models con- 



See [64| for references of low-z techniques for probing the SFR 
density, as well as the non-UV techniques possible at z ~ 2 — 4. 



sidered in the aforementioned paper, the SF2 model at- 
tempts to factor in the uncertainties in the incomplete- 
ness of data sets and the amount of dust extinction at 
early epochs. As such, the SFR density remains roughly 
constant at z > 2. Its form is, 



~dt 



(z) 



0.16 X 



/ exp (3.4z) 



Vexp(3.4z) + 22^ 



(1 + z)3/2 



A/gMpc^'^yr" 



(A2) 



where 



E{z) = y'l7„,o(l + + f^fc,o(l + z)^ + nA{z). (A3) 

Obviously, if studies in the following decade confirm the 
SFR-density trends found in Refs. i65j we would not 
attempt to fit any ET data with the SF2 model. This 
model would need to be updated with a more realistic 
parametrisation. But for now, we adopt the SF2 model 
as a useful ansatz. For the present study, we parametrize 
the SF2 ansatz by making the factors of 3.4 in the nu- 
merator and denominator of Eq. (jA2p variables, labeled 
/?! and respectively. 



Appendix B: A faster calculation of the expected 
detection rate 



The expected detection rate of inspiraling NS-NS bi- 
naries is given by 



Nn = Tx 



' AttD,{zYDh n{z) 
Jo E{z) (l + z) 



X V{M) 



po Dl{z) ( 1.2Mq 



{l + z)M 



5/6' 



dzdM. 
(Bl) 



In our previous study [22!| we found that a simple 
parametrization of the expected detection rate provided a 
good approximation to the slower multi-dimensional inte- 
gration necessitated by Eq. (IBip . However, we now want 
to extend our model-parameter space to a larger number 
of dimensions, for which the simple ansatz method be- 
comes cumbersome. We found in our previous analysis 
that the standard deviation of the NS mass distribution 
had very little impact on the expected detection rate. 
Changing ctns from O.O2M0 to 0.12Mo led to a change 
in the expected detection rate of < 0(1%). 

Although the precision with which we are able to con- 
strain (Tns scales with the number of detections as 1/ a/ZV, 
it is the distribution of detectable systems rather than 
their number that provides information on ctns- If we 
approximate the Gaussian chirp mass distribution by a 
(5-function centered on the mean of the chirp mass distri- 
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bution, we can replace Eq. (jBl[) by a ID integral 



Nd =Tx 



'47TD,{zyDH n(z) 



X Ce 



(1 + z)- 
poi^L(2) f 1.2Mq 



5/6' 



(1 + z)^imJ 



dz. 



(B2) 



This integration can be solved at least an order of 
magnitude faster by standard routines and gives re- 
sults consistent with the full 2D integration procedure. 
We also checked this faster method against the ansatz 
parametrization method, finding that the method used 
in our previous analysis was sufficiently accurate. 
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